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Abstract 

Ubiquitous Van der Waals interactions between atoms and molecules are important for many 
molecular and solid structures. These systems are often studied from first principles using the 
Density Functional Theory (DFT). However, the commonly used DFT functionals fail to capture 
the essence of Van der Waals effects. Many attempts to correct for this problem have been proposed, 
which are not completely satisfactory because they are either very complex and computationally 
expensive or have a basic semiempirical character. We here describe a novel approach, based on 
the use of the Maximally-Localized Wannier functions, that appears to be promising, being simple, 
efficient, accurate, and transferable (charge polarization effects are naturally included). The results 
of test applications are presented. 
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DFT represents a well-established tool to study the structural and electronic properties 
of molecules and condensed matter systems from first principles, and to elucidate complex 
processes such as surface adsorptions, catalytic reactions, and diffusive motions. Although 
current density functionals are able to describe well several systems, at much lower com- 
putational cost compared to other first principles methods, they fail to do so[l[ for the 
description of long-range dispersion effects, generally denoted as Van der Waals (VdW) in- 
teractions, particularly the leading R~ 6 term originated from correlated instantaneous dipole 
fluctuations; the cases where DFT (using, for instance the Local Density Approximation, 
LDA, or the PBeQ functionals) provides reasonable estimates for the interaction energy 
of weakly bound system are actually due to favorable errors or cancellations and should 
therefore be considered accidental. 

In order to overcome this severe deficiency of DFT, two basic strategies have been 
adopted: on one hand, new density functionals or/and relatively complex schemes have 
been proposed that in principle allow for a correct treatment of the VdW interactions 
l|, [^l, I4], Q, ^ Q, y, y], on the other hand several semiempirical approaches 10, 11 1 have 
been developed where an approximately derived R~ 6 term, multiplied by a suitable short- 
range damping function, is explicitly introduced. Although both these approaches have been 
somehow successful, neither of them appears to be entirely satisfactory: in fact, the former 
is generally quite complex and computationally very demanding, compared to a standard 
DFT calculation, while the latter, based on interatomic Cq coefficients (actually dependent 
on the molecular environment of the atoms involved) and empirical fits, turns out to be far 
from generally applicable because it neglects changes in the atomic polarizabilities (which, 
in general, are not additive) and should be tailored to the specific system considered. There- 
fore the development of a practical efficient scheme to include VdW interactions in DFT 
remains a great theoretical challenge. 

In this paper we propose a novel method which allows the efficient calculation of the VdW 
interaction between nonoverlapping fragments, using as input only the ground state electron 
density and the Kohn-Sham (KS) orbitals computed in a conventional DFT approach. 

Crucial to our analysis is the use of the Maximally-Localized Wannier function (MLWF) 
formalism [3], that allows the total electronic density to be partitioned, in a chemically 
transparent and unambiguous way, into individual fragment contributions. The MLWFs 
represent a generalization, for systems characterized by periodic boundary conditions, of 
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the Boys' localized orbitals[l3l| that are commonly used in quantum chemistry; they allow 
for an intuitive interpretation of the bonding properties of condensed-matter systems 



i2h 



and are at the center of the modern theory of polarization [14|. The MLWFs, {w n (r)}, 
are generated by performing a unitary transformation in the subspace of the occupied KS 
orbitals, obtained by a standard DFT calculation, so as to minimize the total spread: 

S = S n = (( w n\r 2 \w n ) - (w n \r\w n ) 2 ) . (1) 



Besides its spread, S n , each MLWF is characterized also by its Wannier-function center 
(WFC); for instance, if periodic boundary conditions are used with a cubic supercell of side 



L, the coordinate x„ of the n-th WFC is defined 



12] as 



x n = -— Imln^le l ^ x \w n ) , (2) 

with similar definitions for y n and z n . If spin degeneracy is exploited, every MLWF corre- 
sponds to 2 paired electrons. Starting from these MLWFs the leading R~ 6 VdW correction 
term can be evaluated using different possible recipes; one of them is described and applied 
in the following. We make the reasonable (at least for insulating systems) assumption JjJ 
of exponential localization of the MLWFs in real space, so that each of them is supposed 
to be an hydrogen-like, normalized, function, centered around its WFC position, r n , with a 
spread S n : 



^(|r-r n |) = -^e-^' r --'. (3) 



3 3 / 4 

Then the binding energy of a system composed of two fragments is given by £& = Eo + Ey^w, 
where Eq is the binding energy obtained from a standard DFT calculation, while the VdW 
correction is assumed to have the form: 

#vdw = - fnl (r n i) -p , (4) 

n,l r nl 

where r n \ is the distance of the n-th WFC, of the first fragment, from the Z-th WFC of the 
second one, the sum is over all the MLFWs of the two fragments, and the C Gn i coefficients 
can be calculated directly from the basic information (center positions and spreads) given 
by the MLFWs. In fact, using for instance the expression proposed by Andersson et al. 
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(see Eq. (10) of ref. Q]) that describes the long-range interaction between two separated 
fragments of matter: 



C 6nl = -J— [ dr [ dr' V /p " (r)P;(r/) = _J_ f dr f dr > ^n(r)w l (r') 

&n 327T 3 / 2 7|r|<r c J\r'\<r' c J p n ( r ) + yj Pl (r') 32tT 3 / 2 7|r|<r c J\r'\<r' c W n (r) + W t (r') 

(5) 



where p n (r) = w^(r) is the electronic density corresponding to the n-th MLWF, C§ n \ is given 
in a.u., and the r c , r' c cutoffs have been introduced 3, J] to properly take into account both 
the limit of separated fragments and of distant disturbances in an electron gas: by equating 
the length scale for density change to the electron gas screening length one obtains: 

6p n (r c ) _ v F [p n {r c )] ^ 



\VPn(r c )\ u p [p n (r c )} ' 
where vp = (37r 2 p n (r)) 1 ^ 3 jm is the local Fermi velocity, and uo p = (4iTe 2 p n (r)/m) 1 ^ 2 is the 
local plasma frequency. By using the analytic form (see Eq. (j3])) of the MLWFs, it is 
straightforward to obtain the cutoff expressed in terms of the MLWF spread: 

r c = S n V?> (0.769 + l/21n(S n )) , (7) 

and to evaluate very efficiently the multimensional integral of Eq. (jSJ). For instance, in 
the test case of 2, distant H atoms, using the well known (unperturbed) analytic H atom 
wavefunction, the above formula gives C§ = 6.41 a.u. to be compared to the reference 
literature value of 6.50 a.u. 

In Eq. (jHJ), if the electronic density corresponding to every MLWF is multipled by 2, the 
C§ni coefficients increase by a y/2 factor; therefore it appears reasonable to assume that, 
when each MLWF describes 2 paired electrons (spin degeneracy), C 6n ; has to be multipled 
by y/2. This is also supported by the fact that, in the Slater-Kirkwood approximation for 
estimating the Cq coefficients, the effective number of electrons is smaller than the number 
of valence electrons, and it is 1.42 ~ y/2 in the case of the He atomjlf]], whose DFT ground 
state is just given by 2 paired electrons in the lowest-energy KS orbital. 

In Eq. (J4]) f n i{r) is a damping function which serves to cutoff the unreasonable behavior 
of the asymptotic VdW correction at small fragment separations. For it we have chosen a 



form 
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17l | with parameters directly related to the MLWF spreads: 
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fm{r) 



1 

1 + exp(—a(r / R s — 1)) 



(8) 



where [171] a ~ 20 (the results are almost independent on the particular value of this param- 
eter), and R s = -Rvaw + -^vdw * s the sum °f the VdW radii of the MLWFs, which, following 
Grimme et al. JJTj], are determined as the radii of the 0.01 electron density contour; using 



Eq. (j3J) one easily obtains that: 



R YdW = (1-475 - 0.8661n(S n ))S n . (9) 

The damping function effectively reduces the VdW correction to zero typically below 2 A; 
at intermediate distances a minimum in the VdW potential exists that usually lies slightly 
below the sum of the corresponding VdW radii, R s . Note that the above recipe resembles 



that proposed in ref. 



181 ]. where the long-range electron-electron interaction is separated by 



the short-range one, using a single parameter describing the physical dimensions of a valence 
electron pair. 

The Eq binding energy can be obtained from a standard DFT calculation (we have used 
the CPMd[i9] and I/-ESPRESSO 20] ab initio packages), using t re Generalized Gradient 
Approximation (GGA) in the revPBE flavor \^\. This choiceja, 1t| is motivated by the 
fact that revPBE is fitted to the exact Hartree-Fock exchange, so that the VdW binding, a 
correlation effect, only comes from the VdW correction term, as described above, without 
any double-counting effect (for instance, LDA, or other GGA functionals, such as PBE, 
predict substantial binding in rare gas dimers, due to a severe overestimate of the long-range 
part of the exchange contribution ]8]). The evaluation of the VdW correction as a post- 
standard DFT perturbation, using the revPBE electronic density distribution, represents an 
approximation because, in principle a full self-consistent calculations should be performed; 
however recent investigations 22| on different systems have shown that the effects due to the 
lack of self-consistency are negligible (this is reasonable because one does not expect that 
the rather weak and diffuse VdW interaction substantially changes the electronic charge 
distribution). 

The VdW correction scheme described above can be refined by considering the effects due 
to the anisotropy of the MLWFs, and distinguishing between contributions along (or orthog- 
onal to) the fragment-fragment direction (details will be published elsewhere [231] ). Moreover, 
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also higher-order term VdW corrections, involving the Cg, C 10 ,-~ coefficients, could be eas- 
ily included. Clearly, in the present method, the evaluation of the VdW corrections to the 
interfragment forces is trivial, thus allowing an easy implementation in standard geometry 
optimization calculations or Molecular Dynamics simulations. Remarkably, the whole pro- 
cedure of generating the MLWFs and evaluating the VdW corrections represents a negligible 
additional computational cost, compared to that of a standard DFT calculation. 

We have applied the new method to selected dimers among typical VdW-bonded systems: 
Ar 2 , N 2 -N 2 ("T-shaped"), CH 4 -CH 4 , C 6 H 6 -C 6 H 6 ("sandwich-shaped"), C 6 H 6 -Ar, C0 2 -C0 2 , 
and also a mixed (H-bonded/VdW-bonded) complex, C 6 H 6 -H 2 0. 

In the Tables I-III we report our computed binding energy (a positive value indicates 
unbound complexes), equilibrium characteristic interdimer distance, and fragment-fragment 
effective C6 coefficient, which is defined as J2n,i^eni (in fact, neglecting the differences in 
the spreads of the MLWFs and in the interfragment distances between the WFCs, -Eydw ~ 
— / (iLnjCeni) /R 6 , where R is the interdimer distance). These values are compared to the 
most reliable (to our knowledge) reference literature corresponding data, which are often 
spread over a relatively large range (comparison with literature C6 coefficients should be 
taken as purely indicative because of different assumed definition). As can be seen, the 
general performance of the method is quite satisfactory; in fact, the improvement achieved 
by including the VdW correction, with respect to the pure revPBE results, is dramatic, 
even in the case of a mixed complex, such as CeH 6 -H 2 0, where some fraction of the binding 
energy is already given by the standard DFT calculation. In Fig. 1 we show the effect of 
the inclusion of the VdW correction on the behavior of the binding energy of the Ar 2 dimer, 
plotted as a function of the Ar-Ar distance, and compared to the reference equilibrium value; 
in Ar the 8 valence electrons of each atom are described by 4 MLWFs (spin degeneracy is 
exploited), whose WFCs are tetrahedrally located around the Ar ion. 

In the case of the equilibrium characteristic distances, the more substantial deviation 
from the reference results can easily be explained by the fact that the potential energy 
curves for weakly-bonded systems are typically very shallow. Inspection of Table I shows 
that anisotropy effects do not much affect the binding energy estimates, but for the case of 
the "sandwich-shaped" benzene dimer, where anisotropy correction leads to a value much 
closer to the reference data: this behavior comes as no surprise, since the planar geometry of 
the two fragments clearly induces strong anisotropy in the computed MLWFs. In Table I, by 
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considering the MLWF anisotropy, the binding energies are always decreased and lower than 
the reference values; this behavior is probably due to the neglect of higher-order contributions 
to the VdW correction, such as the —C$/R s term (dipole-quadrupole interaction), which 



24]. 



should be included to have a very accurate estimate of the binding energy 

We have also applied our technique to the case of graphite, where the optimal interlayer 
distance, found to be 8.69 A with the standard DFT-revPBE approach, becomes 6.33 A by 
including VdW effects, that is much closer to the experimental value (6.70 A). 

In conclusion, we have presented and applied in test cases a technique suitable to describe 
VdW effects in the framework of standard DFT calculations. The technique is based on the 
generation of the MLWFs and naturally describes changes in the electronic density distri- 
butions of the fragments due to the environment, for instance related to charge polarization 
effects: in fact these changes are easily described in terms of changes in the location of the 
centers and in the spreads of the MLWFs. The results of the method, which is simple to 
be implemented and not expensive computationally, are quite satisfactory and promising, 
also considering that a large area for future improvements exists: in fact, different, more 
sophisticated schemes to utilize the MLWFs could be developed and/or improved reference 
DFT functionals, with respect to revPBE, could be adopted. 

We thank F. Ancilotto, M. Boero, F. Toigo, and F. Valencia for useful discussions. We ac- 
knowledge funding from Padua University project n. CPDA033545 and PRIN 2004 ("NAN- 
OTRIBOLOGIA" project), and allocation of computer resources from INFM "Progetto 
Calcolo Parallelo". 
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TABLE I: Binding energy, in meV, computed using the standard DFT-revPBE calculation, Eq and 
including the VdW correction, Eq + £VdW> compared to available literature data; in parenthesis 
values, computed taking anisotropy effects into account, are reported. 



system 


Eq 


Eq + -EVdW 


ref 


Ar-Ar 


-1.7 


-11.9 (-9.5) 


-12.3 


N 2 -N 2 


-2.9 


-11.1 (-10.8) 


-13.3 


CH4-CH4 


-2.1 


-11.7 (-9.9) 


-23<->-14 




+7.1 


-252.7 (-144.1) 


-142^-61 


C 6 H 6 -Ar 


-2.4 


-65.7 (-53.5) 


-65 


C0 2 -C0 2 


-16.2 


-54.0 (-47.9) 


-69^-59 


CeH6-H 2 


-40.4 


-131.2 (-121.6) 


-169^-137 



TABLE II: Equilibrium characteristic interdimer distance, in A, computed using the standard DFT- 
revPBE calculation, Rq and including the VdW correction, R\dw, compared to available literature 
data; in parenthesis values, computed taking anisotropy effects into account, are reported. 



system 


Rq 


i?VdW 


ref 


Ar-Ar 


4.67 


4.03 (4.07) 


3.76 


N 2 -N 2 


5.05 


4.37 (4.37) 


4.03 


CH4-CH4 


4.70 


4.23 (4.25) 


3.60^4.27 






3.45 (3.45) 


3.80^3.90 


C 6 H 6 -Ar 


4.79 


3.57 (3.57) 


3.41 


C0 2 -C0 2 


3.86 


3.49 (3.49) 


3.60 


CeH6-H 2 


4.23 


3.17 (3.17) 


3.40^3.50 



TABLE III: Fragment-fragment effective Cq coefficient (see text for the definition), in a.u., com- 
puted using the standard DFT-revPBE with the VdW correction, compared to available literature 
data; in parenthesis values, computed taking anisotropy effects into account, are reported. 



system 


c 6 


ref 


Ar-Ar 


92.5 (74.7) 


64.3^65.5 


N 2 -N 2 


90.3 (95.6) 


73.4 


CH4-CH4 


103.0 (101.0) 


118.0^130.0 


C6H6-C6H6 


2930.0 (2460.0) 


1723.0 


C 6 H 6 -Ar 


490.0 (448.0) 


330.1 


C0 2 -C0 2 


187.0 (172.0) 




CeH 6 -H 2 


323.0 (299.0) 


208.5^277.4 




(A 9m ) X§J9U9 Suxpuxq 



FIG. 1: Binding energy of the Ar2 dimer, as a function of the Ar-Ar distance, using the standard 
DFT-revPBE calculation and including the VdW correction. 



